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ABSTRACT 

Heat  transfer  rates  for  laminar,  convective  heat  transfer  in  the 
entrance  region  between  parallel  plates  were  investigated.   The  hydro- 
dynamic  solution  due  to  Bodia  f"2l  was  used  in  the  solution  of  the  energy 
equation  in  finite  difference  form  on  a  digital  computer.   The  thermal 
boundary  conditions  include;  constant  heat  inputs  constant  wall  tempera- 
ture, one  wall  constant  temperature  and  one  wall  insulated,  and  constant 
but  different  wall  temperatures  on  the  upper  and  lower  walls. 

The  approximate  integral  methods  of  Siegel  and  Sparrow  T5l,  TT]   are 
shown  to  produce  results  that  are  in  close  agreement  with  the  solutions 
in  the  present  analysis  for  the  constant  heat  input  and  constant  wall 
temperature  cases. 

The  scope  of  the  finite  difference  solution  is  limited  to  a  narrow 
range  of  Prandtl  numbers  near  unity,  due  to  the  small  grid  sizes  required 
for  convergence. 
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I.  INTRODUCTION 

Efficient  design  of  compact  heat  exchangers  such  as  are  employed  for 
gas  turbine  regenerators,  among  other  applications,  requires  maximum  utili- 
zation  of  the  high  transfer  rates  available  in  the  entrance  region.   In 
this  region  the  velocity  and  temperature  profiles  of  the  entering  fluid  are 
undergoing  rapid  transition  from  their  uniform  distribution  at  entrance  to 
the  fully  established  distributions  encountered  farther  downstream.   The 
well  known  Graetz  solution,  in  which  only  the  temperature  development  is 
considered,  does  not  yield  adequate  results  when  the  velocity  and  tempera- 
ture distributions  are  simultaneously  developing.   This  is  particularly 
true  if  the  entrance  region  represents  a  considerable  percentage  of  over- 
all length  and  the  two  profiles  develop  at  approximately  the  same  rates. 
This  situation  is  frequently  encountered  with  gases  where  the  Prandtl 
number  is  approximately  one. 

The  purpose  of  the  present  analysis  is  to  present  a  finite  difference 
solution  for  the  heat  transfer  rates  between  parallel  flat  plates  for 
several  different  boundary  conditions  corresponding  to  various  constant 
wall  temperatures  and  constant  heat  input  configurations.   The  approach 
to  the  problem  was  suggested  by  Miller  Til,  which  is  adapted  from  the  num= 
erical  procedures  applied  to  the  entrance  region  of  the  circular  tubes  by 
Kays  T2"). 

While  an  analysis  such  as  the  present  one  can  hardly  compete  for 
simplicity  with  approximate  treatments  such  as  the  Graetz  analysis,  or 
the  several  energy  integral  analyses  commonly  employed  in  engineering  design, 
it  does  provide  an  exact  solution  for  a  limited  number  of  cases  which  serves 
as  a  standard  with  which  the  accuracy  of  the  various  approximate  methods 
may  be  assessed. 
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II.  ANALYSIS 

A.   Governing  Equations 

The  governing  equations  are  the  energy,  momentum  (Navier-Stokes)  and 
continuity  equations.   In  order  to  reduce  the  complexity  and  coupling  of 
the  equations,  the  following  assumptions  are  made. 
The  flow  is  assumed  to  be: 

1.  steady s  two-dimensional 

2.  laminar 

3.  incompressible 

In  addition  it  is  assumed  that: 

1.  thermal  diffusivity  (ex)  is  constant 

2.  convective  heat  transfer  is  large  compared  to  radiation,  axial 
conduction,  and  viscous  dissipation. 

The  governing  equations  in  reduced  form  may  then  be  expressed  as: 

Energy: 

St  ,   at     _    d2t 

u  —  +  v  -r—  =  oc  «-  n\ 

dx     dy       dy2  UJ 


Momentum: 


du  _l   ^u      1  dp  ,   ^2 

u  — —  +  v  t—  •  -  —  -2£  +  u  2  u 

SX       By         P   oX       £yT 


(2) 


Continuity: 


|H..+  |v  ,0  (3) 

dx   dy 

Due  to  the  assumption  of  incompressibility  the  above  equations  are 
no  longer  coupled.   The  temperature  development  remains  dependent  upon 
the  velocity  profile,  although  dependence  of  velocity  on  temperature  has 
been  removed. 


The  boundary  conditions  associated  with  equations  (1),  (2),  and  (3) 
for  flow  between  parallel  flat  plates  with  uniform  velocity  and  temperature 
at  entrance  are: 

1.  At  the  entrance  (x  ■  0) 
t--  to  -Ci 

u  «  uQ  -  C2 

v  ■  v0  =  0 

2.  At  the  walls  (y  ■  +  d) 
v  =  vw  =  0 

u  =  u  =0 
w 

fc  ■  fcw 


3.  On  centerline 
v  =  vc  -  0 

4.  Applied  thermal  boundary  conditions 

(1)  symmetric  cases  (with  respect  to  duct  centerline) 

a.  constant  wall  temperature 

t  »  C~ 
w    3 

b.  constant  heat  input 

(2)  asymmetric  cases 

a.  one  wall  constant  temperature,  one  wall  insulated 

Cwi  =  C5 

b.  Constant  wall  temperature 

Cwi  =  c6 


.    . 


.    . 


tw2    C7 
where:   C^  +   C7 

Solution  of  the  energy  equation  requires  that  a  hydro dynamic  solution  be 
first  obtained. 
B.   Hydrodynamic  Solution 

Several  hydrodynamic  solutions  for  the  entrance  region  are  available. 
Solutions  have  been  presented  by  Bodoia  T3lj  Schlichting  T4]s  Schiller  r5~ls 
and  Han  ^61.   The  approximate  solutions  of  Schiller  and  Han  have  the  advantage 
of  being  analytic.,  and  therefore,  less  cumbersome  to  use.   They  are,  however, 
inherently  less  accurate  than  the  numerical  solutions  presented  by  Bodoia  and 
Schlichting. 

Bodia's  solution  was  obtained  by  evaluation  of  the  Prandtl  boundary  layer 
equations  in  finite=dif ference  form.   Velocity  profiles  were  obtained  at  given 
streamwise  locations  by  utilizing  matrix  methods  to  simultaneously  solve  a 
column  array  of  momentum  expressions  normal  to  the  flow. 

A  comparison  of  Schlichting' s  series  solution  and  the  Bodoia  solution 
indicates  thats  in  addition  to  the  velocity  gradient  discontinuity  in  the 
Schlichting  solution  where  the  upstream  and  downstream  solutions  are  joined, 
Schlichting' s  representation  of  the  pressure  gradient,  based  only  on  the 
centerline  velocity,  results  in  a  core  velocity  gradient  which  is  too  large 
and  a  wall  velocity  gradient  which  is  too  small  f3l. 

Comparisons  of  the  dimensionless  pressure  decrement  versus  length  due 
to  Bodoia,  Schlichting  and  Han  are  presented  in  Fig.  2.   Velocity  distribution 
comparisons  between  Bodoia  and  Schlichting,,  and  Schlichting  and  Schiller,  are 
shown  in  Figs.  3  and  4. 

The  Bodoia  solution  was  adopted  for  the  present  analysis. 
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C.   Energy  Solution 

Equations  (1)  and  (3)  may  be  written  in  dimensionless  form  by  introduc- 
ing the  following  dimensionless  variables! 

T  -  t/t  X  -  22^r~     -  x/d  U  ■  u/uQ 

d  uo     Re 
Y  -  y/d  V  =  v/uQ 

In  dimensionless  form,  the  equations  are  then: 
Energy: 

oX      9Y     PR  dY7  k  ; 

Continuity: 

The  following  finite  difference  approximations  are  then  introduced: 

3T,^3  Tx+lay  -  Tx,y  „. 

dX   Ax       A  X 

ST  =  AT  =  Tx,y+1  -  Tx,y-1 

SY   AY   2~A~y (7> 


2 

g  T  _  Tx,y+1  -  2Tx,y  +  Tx,y°l 
—7-  - 


(8) 
SX  A  Y2 

SU  =  Aq  _  Ux+l,y-Ux-l,y  m 

SX   AX      2  A  X 

dV  _  Ay  ..  Vx,y~Vx,y-l 

oY   AY   Ay (10) 

Note  that  (6)  is  a  forward  difference  and  (10)  is  a  backward  difference, 
whereas  the  remaining  equations  are  of  the  central  defference  type.   These 
have  been  expressed  in  this  manner  for  simplicity  and  convenience.   The  larger 
the  truncation  error  normally  associated  with  the  forward  or  backward  difference 
is  not  significant  here  since  practical  grid  dimensions  require  that: 
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AY  »  AX 

and,  due  to  the  nature  of  the  flow? 

Au    Av 
Ax  >>  AY 

Introducing  (6)  through  (10)  into  (4)  and  (5)  and  solving  for  the  ap= 
propriate  term  yields: 


Tv+1  _,  _  I  AX     .  AXReV  ,  m    , 
Tx+1'y   L  UPRAY2  +  ^2AY"  J  TX'y_1  + 


i  .    2AX 

L  l  "  upray7  J  Tx>y+1  + 


I    AX     AXReV   |  m    . ,  /11X 

L  UPRAY7  "  U2AY~  J  TX'y+1  (U) 


Ay 

Vx,y  =  VXjy-l  +RelAl  [  Ux_1'y  "  Ux+1'y  J        (12) 
which  contribute  the  computing  equations. 

Note  that  the  Reynolds  number  (Re) s  although  appearing  in  the  above 
equations,  remains  only  for  convenience  and  could  be  eliminated.   It  is  no 
longer  a  parameter  in  the  set  of  equations  haveing  been  included  in  the 
dimensionless  length  variable  (X). 

A  mixed  mean  temperature  (t  )  may  be  expressed  in  the  following  way: 


(o.    putdA 


"m   Jc    pudA  (13) 

Removing  the  constants  from  within  the  integral ,  and  simplifying,  a 
dimensionless  mixed  mean  temperature  may  be  written  in  finite  difference 
form: 


y-1 
m     L    Ux,y  Tx,y 


m   y=o 

(14) 


y-1 

L  Ux>y 

y-o 

Due  to  the  no-slip  condition  imposed  upon  viscous  flows,  heat  trans- 
fer to  or  from  the  fluid  at  the  wall  may  be  expressed  in  terms  of  con- 
duction through  a  thin  film  of  stagnant  fluid.   Thus: 

<•"  -  </A  "  *  (f  ^  (15) 

which  must  equal  the  heat  flux  due  to  convection, 

q"  «  q/A  -  h(tw  -  tm)  (16) 

Defining  the  local  Nusselt.  number  as: 

N   =  M  (17) 

k 

Combining    (14)    and   (15), 

N„„  =     V5yA-o  (18) 


ux 


fcv-   Cm 


A  Nusselt  number  based  on  the  hydraulic  diameter  (d)  instead  of  half 
the  gap  width  (d)  in  finite  difference  form  becomes: 

4(t  -  T 


\   w    w-Ay) 


NUx  -     T  J* <U> 

w    m 


The  hydraulic  diameter  was  introduced  in  (19)  to  enable  comparisons 
with  previous  solutions.   The  conversion  is: 

D  =  4d 

In  addition,  for  the  asymmetric  constant  wall  temperature  case,  it 
is  convenient  to  use  a  dimension less  parameter  other  than  the  Nusselt 


number.   A  dimensionless  heat  flux  parameter  may  be  defined  by  using  (15) 


above: 

kt. 


q*  = 


fe-^^-if^J-^ 


D.  Thermal  Boundary  Conditions  in  Finite  Difference  Form  Constant  Wall 
Temperature  (Symmetric  Case).   Since  the  temperature  distribution 

is  symmetrical  with  respect  to  the  centerline,  computations  were  performed 
in  only  the  upper  half  of  the  duct.   A  coordinate  system  was  adapted  with 
the  X  axis  on  the  centerline  of  the  channel.   The  thermal  boundary  condi- 
tion at  the  wall  was  applied  as: 

t  /t  =  T  =  N 
wo    w 

where  N  is  a  positive  integer. 

Constant  Heat  Input  (Symmetric  Case).   With  the  axis  system  as  above, 
the  slope  of  the  thermal  profile  was  maintained  constant  at  the  wall  by 
adding  a  constant  to  the  calculated  temperature  one  step  from  the  wall, 
thus: 

T  =  T  =A  +  M 

w    w  Ay 

where : 

m   =  (At /ay)  Ay 

w 

Constant  Wall  Temperature  with  One  Wall  Insulated.   For  this  case  the 
coordinate  system  was  redefined  so  the  X  axis  coincided  with  the  lower  wall, 
The  lower  wall  boundary  condition  was  taken  to  be: 

t   /t  »  T   »  N 

Wl   O      Wi 

while  zero  slope  at  the  upper  wall  was  obtained  by  equating  the  wall 
temperature  to  the  calculated  temperature  one  step  from  the  wall: 


V'o  =  Tw2  =  Tw2  =  AY 

Constant  Wall  Temperature  (Asymmetric  Case).   With  the  axis  system 
the  same  as  for  the  insulated  wall  case,  thermal  boundary  conditions  were 
established  as: 

Vfco  =  T"l  =  Nl 

tW2' to  =  ^w2  =  ^2 
where  N-^  and  N2  are  positive  integers  and  N^  N2^ 
Noting  that: 


trT-l    -   t 
wi          o 

= 

^l^o   ■    l 

Nl  - 1 

tW2   '     to 

V*o  "   L 

N2   -    1 

Values  of  N,  -  1/N?  -  1  of  1/2,  1/3,  1/4  were  investigated. 
E.   Computational  Procedure 

The  finite  difference  solution  was  obtained  by  using  a  CDC  1604 
digital  computer  to  solve  the  equations  at  intervals  of  AX  ■  0.001  and 
AY  ■  0.1.   Velocity  distributions  were  obtained  by  interpolating  plot= 
ted  curves  of  Bodoia's  tabulated  results. 

Since  the  entire  solution  is  dependent  upon  the  first  calculations 
at  the  entrance,  starting  values  were  obtained  by  reducing  the  grid  size 
to  Ax  ■  0.0001  and  AY  ■  0.05  for  the  first  20  streamwise  stations  from 
X  -  0  to  0.0020. 

The  continuity  equation  (2)  was  evaluated  and  the  cross-wise  velocity 
component  (V)  was  introduced  into  the  energy  equation  (11)  at  each  grid 
point.   The  energy  solution  was  then  obtained  at  intervals  of  AY  between 
centerline  and  the  wall  at  each  streamwise  station. 


Solutions  for  each  thermal  boundary  condition  were  obtained  for 
Prandtl  numbers  of:   0.5,  0.7,  1.0,  1,6,  3.2  and  10. 

Solutions  for  Prandtl  numbers  less  than  0.5  are  excluded  by  the  grid 
size.   This  may  be  seen  from  inspection  of  the  coefficient  of  the  second 
temperature  term  in  (11),  which  indicates  the  interdependence  of  PR,  AX, 
and  AY.   In  order  to  obtain  convergence,  PR  must  be  large  enough  to  pre- 
vent this  term  from  becoming  negative.   Thus  low  Prandtl  numbers,  ap- 
proaching those  of  liquid  metals,  cannot  be  effectively  handled  by  finite- 
difference  methods  owing  to  the  large  number  of  calculations  required  by 
the  small  grid. 

The  computer  program,  written  in  Fortran  60,  is  included  in  the 
Appendix. 
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Ill  RESULTS 

Tables  I  through  VI  contain  computed  and  extrapolated  values  of  the 
heat  transfer  rates  as  a  function  of  distance  from  the  entrance  for  each 
of  the  thermal  boundary  conditions  for  a  Prandtl  number  of  0.7.  Figures 
5  through  29  contain  curves  of  the  transfer  rate  versus  downstream  posi= 
tion  and  temperature  and  velocity  profiles  for  Prandtl  numbers  of  0.5, 
0.7,  1.0,  1.6,  3.2  and  10. 

With  the  exception  of  the  Asymmetric  Constant  Wall  Temperature  Case, 
which  is  expressed  in  terms  of  the  dimensionless  heat  flux  parameter,  q*; 
the  heat  transfer  rates  are  presented  as  local  Nusselt  number  variation 
with  the  length  parameter,  L.   Whereas,  the  computations  were  carried  out 
using  the  dimensionless  length  parameter  X,  which  is  normally  employed 
for  hydrodynamic  developments;  L,  based  on  the  hydraulic  diameter  is 
customarily  employed  in  heat  transfer.   The  conversion  is: 

X  -  16L 

The  use  of  all  finite  starting  length  in  the  finite  difference  equa- 
tions results  in  a  finite  initial  value  for  the  Nusselt  number  at  X  =  0 
rather  than  the  theoretically  infinite  value  at  the  entrance.   This  causes 
a  nearly  flat  slope  near  the  entrance  and  an  artificial  inflection  point 
in  the  Nusselt  number  versus  length  curve.   The  inflection  point  occurs 
progressively  further  downstream  with  the  increasing  Prandtl  number, 
thereby  limiting  the  useful  span  of  the  curves.   A  reduction  of  grid  size 
from  AX  =  .001,  AY  =  . 1  to  ^X  =  .0001,  AY  -  .05  was  observed  to  double 
the  initial  Nusselt  number  from  40  to  80,  and  move  the  inflection  point 
upstream.   This  increased  the  useful  span.   However,  the  approach  obviously 
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reaches  a  practical  limit  very  rapidly  due  to  the  number  of  calculations 
and  the  accuracy  required  of  the  hydrodynamic  solution  very  close  to  the 
entrance. 

In  general s  for  the  three  thermal  boundary  conditions  in  which  the 
Nusselt  number  was  used  as  a  heat  transfer  parameter,  the  curves  for 

Prandtl  numbers  of  0.5,  0.7,  1.0  and  1.6  behave  as  expected  for  values  of 

-A.  -3 

L  >10~  .   For  Prandtl  numbers  of  3.2  and  10,  L  ■  10   is  the  approximate 

lower  limit  of  validity. 

Constant  Heat  Input.   A  comparison  of  local  Nusselt  numbers  obtained 
in  the  present  analysis  with  those  reported  by  Siegel  and  Sparrow  [*7l,  and 
Han  T6],  contained  in  Fig.  5,  indicate  close  agreement  with  the  results  of 
Siegel  and  Sparrow.   Their  solution  employed  a  simplified  energy  equation 
in  which  the  temperature  distribution  in  the  boundary  layer  was  expressed 
as  a  series  of  polynomials  in  the  transverse  coordinate  using  the  down- 
stream station  as  a  parameter.   The  hydrodynamic  solution  employed 
Schiller's  approximation. 

The  comparison  with  Han's  integro-numerical  solution  is  not  as  close 
as  might  appear  in  Fig.  5  since  this  curve  is  for  a  Prandtl  number  of  0.8 
and  should  lie  above  the  other  two  curves  which  are  for  a  Prandtl  number 
of  0.7. 

The  computed  curves  of  local  Nusselt  number  versus  length  contained 

in  Fig.  6,  for  Prandtl  numbers  of  0.5,  0.7,  1.0,  1.6,  3.2  and  10,  may 

=  3 
be  approximated  within  10  percent,  for  0.5  «£PR£l.6,  and  L  >10   by  the 


following: 


(PR/L) 

N   =  8.24  +  iP-I**6. r 

ux  1  +  .0178  (PR/L)- 6 
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• 


Velocity  and  temperature  profiles  for  this  boundary  condition  are  con- 
tained in  Figs.  13  through  16. 

One  Wall  Constant  Temperature  °  One  Wall  Insulated.   The  results  of  this 

computation  are  presented  in  Fig.  17.   The  local  Nusselt  number  may  be  approx- 

=  3 
imated  within  10  percent  for  0.5^PR^1.6,  and  L  >10   by: 


PR 

,AR,       .0155rK/L) 
Nu   ■  4.84  +  Sp ' 

x  1  +  .012(^K/L)° 


The  velocity  and  temperature  distributions  are  illustrated  in  Figs. 
18  through  21. 

Constant  Wall  Temperature  (Asymmetric  Case).   Dimensionless  heat  flux 

values  for  this  boundary  condition  are  presented  in  Figs.  22,  23  and  24 

for  Nl"1  ratios  of  1/2,  1/3,  and  1/4. 
N2-l 

Note  that  the  upper  (hotter)  wall  curve  is  smoothly  asymptotic  to  the 

value  predicted  for  the  steady  state  condition.   The  lower  curve  exhibits  a 

mush  shallower  initial  slope,  passes  through  an  inflection  pointy  and  thena 

approaches  a  negative  asymptotic  value  of  the  same  magnitude  and  at  the 

same  rate  as  the  upper  curve. 

Velocity  and  temperature  profiles  for  N,-l        are  presented  in 

~-r=   1/4 
Figs.  25  through  28.  2 

The  behavior  of  the  local  Nusselt  number  for  this  asymmetrical  case 

is  depicted  in  Fig.  29.   The  lower  curve 3  representing  the  local  Nusselt 

number  variation  for  the  lower  (cooler)  wall,  does  not  depict  the  actual 

physical  situation. 
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Both  the  denominator  and  the  numerator  of; 


4(T   -  T  A  ) 


Ay 

Nux H"T <18> 

x        w    m 

approach  zero  at  different  rates.   The  curve  tends  to  positive  infinity 
as  the  denominator  approaches  zero,  returns  from  negative  infinity,  and 
becomes  zero  as  the  numerator  goes  to  zero;  and  finally  approaches  the 
actual  steady  asymoptotic  value. 
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IV.   CONCLUSIONS 
From  the  results  of  this  study  it  may  be  concluded  that: 

1.  The  finite  difference  method,  as  employed  in  this  analysis,  yields  re- 
sults which  are  consistent  with  previous  solutions  for  other  geometries  [2] 
and  which  are  asymptotic  to  the  known  fully  developed  values.   The 

method  is,  however,  limited  to  a  small  range  of  Prandtl  numbers  near 
unity  by  the  small  grid  sizes  necessary  outside  this  range.   For  small 
Prandtl  numbers  the  grid  must  be  very  small  to  obtain  convergence, 
while  for  large  Prandtl  numbers  practical  grid  dimensions  result  in 
significant  inaccuracies  near  the  entrance  due  to  the  "finite  starting 
length". 

2.  The  results  of  the  approximate  integral  methods  suggested  by  Siegel  and 
Sparrow  T7]  for  constant  heat  input,  and  by  Sparrow  f5l  for  constant 
wall  temperature,  compare  very  favorably  with  the  finite  difference 
solutions  reported  here.   This  additional  substantiation  of  the  methods 
suggested  by  these  authors  is  significant  because  of  the  practical 
utility  of  their  approach,  which  requires  relatively  few  calculations, 
is  capable  of  handling  a  wide  range  of  Prandtl  numbers,  and  can  pre- 
dict heat  transfer  rates  very  close  to  the  entrance. 
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TABLE  I 
CONSTANT  HEAT  INPUT 
PR  =  0.7 

ux  w 

.0016  .0001  46.2  1.09 

.0032  .0002  30.0  1.14 

.0064  .0004  23.7  1.18 

.0128  .0008  18.0  1.24 

.0160  .0010  16.6  1.26 

.0320  .0020  12.9  1.35 

.0640  .0040  10.3  1.47 

.1280  .0080  8.87  1.62 

.1600  .0100  8.62  1.67 
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TABLE  II 
CONSTANT  WALL  TEMPERATURE 
(SYMMETRIC  CASE) 
PR  =  0.7 

X  L 

.0016  .0001 

.0032  .0002 

.0064  .0004 

.0128  .0008 

.0160  .0010 

.0320  .0020 

.0640  .0040 

.1280  .0080 

.1600  .0010 


Asymptotic  to 


% 

34. 

,0 

26, 

,0 

19, 

,9 

15, 

,3 

13, 

,7 

10, 

,5 

8, 

,58 

7, 

,78 

7, 

,63 

7, 

,60 
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TABLE  III 
ONE  WALL  CONSTANT  TEMPERATURE  -  ONE  WALL  INSULATED 

PR  =  0.7 

X                     L  Nux  Tw 

.0016                  .0001  33.0  1.000 

.0032                  .0002  24.0  1.000 

.0064                  .0004  19.2  1.000 

.0128                  .0008  13.0  1.000 

.0160                  .0010  11.6  1.000 

.0320                  .0020                9.05  1.000 

.0640                  .0040                7.20  1.000 

.1280                  .0080  5.93  1.001 

.1600                  .0100                5.55  1.021 

.4000                  .0250                4.96  1.139 

Asvnrototic   to      4.84  2.000 
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X 

.0032 

.0064 

.0128 

.0160 

.0320 

.0460 

.1280 

.1600 

.3200 

.6400 

.9600 


TABLE  IV 

CONSTANT  WALL  TEMPERATURE 

(ASYMMETRIC  CASE) 

Nj  -  1  =  1/2    PR  =  0.7 
No  -  1 


L 

* 

91 

.0002 

26.2 

.0004 

17.3 

.0008 

11.9 

.0010 

10.7 

.0020 

8.08 

.0040 

6.18 

.0080 

4.73 

.0100 

4.24 

.0200 

2.24 

.0400 

-.19 

.0600 

-1.23 

Asymptotic  to 

-2.00 

* 

q2 

57. 

,5 

41. 

4 

28. 

,7 

25, 

,9 

18, 

,7 

13, 

,7 

10, 

,1 

9, 

,15 

6, 

,35 

3, 

,81 

2, 

,76 

+2, 

,00 

20 


X 

.0032 

.0064 

.0128 

.0160 

.0320 

.0640 

.1280 

.1600 

.3200 

.6400 

.9600 


TABLE  V 
CONSTANT  WALL  TEMPERATURE 

(ASYMMETRIC  CASE) 
NX  -  1  =  1/3     PR  »  0.7 


N2  -  1 

L 

* 

11 

.0002 

28.0 

.0004 

17.3 

.0008 

11.8 

.0010 

10.8 

.0020 

8.09 

.0040 

6.18 

.0080 

4.69 

.0100 

4.15 

.0200 

1.65 

.0400 

-1.58 

.0600 

-2.97 

Asymptotic  to 

-4.00 

* 

q2 

86. 

,0 

62, 

,0 

43, 

,0 

38, 

,9 

28, 

,1 

20, 

,6 

15, 

,2 

13, 

,8 

9, 

,84 

6, 

,41 

5. 

,02 

+4. 

00 
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X 

.0032 

.0064 

.0128 

.0160 

.0320 

.0640 

.1280 

.1600 

.3200 

.6400 

.9600 


TABLE  VI 
CONSTANT  WALL  TEMPERATURE 

(ASYMMETRIC  CASE) 
Nx  -  1  -  1/4     PR  =  0.7 


N2  -  1 

L 

* 
11 

.0002 

28.2 

.0004 

16.8 

.0008 

12.2 

.0010 

10.5 

.0020 

8.09 

.0040 

6.17 

.0080 

4.66 

.0100 

4.06 

.0200 

1.06 

.0400 

-2.96 

.0600 

-4.71 

Asymptotic  to 

-6.00 

q2 

16. 

,0 

82, 

,9 

57. 

,7 

50. 

,4 

37, 

,5 

27, 

,5 

20. 

,3 

18, 

,4 

13, 

,3 

9. 

,03 

7, 

,27 

+6. 

,00 

22 
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TEMPERATURE     AND  VELOCITY     PROFILES 

CONSTANT    HEAT    INPUT 
i  X  =  0.005  L=0.0003I25 
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TEMPERATURE     AND    VELOCITY     PROFILES 
X=0.Q50  L  =  0.003125 
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TEMPERATURE    AND    VELOCITY      PROFILES 
X=  0.250  L  =0.0156 
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I  CONSTANT     WALL     TEMPERATURE 
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FIGURE    18 


TEMPERATURE    AND   VELOCITY    PROFILES 
ONE   WALL    CONSTANT   TEMPERATURE 
ONE    WALL     INSULATED 
X=  0.005  L  =  0.0003125 
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FIGURE      19 

TEMPERATURE    AND    VELOCITY     PROFILES 
ONE    WALL  CONSTANT    TEMPERATURE 
ONE     WALL     INSULATED 
X=0,050  L  =0.003125 
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FIGURE      20 


TEMPERATURE     AND    VELOCITY      PROFILES 
ONE    WALL    CONSTANT    TEMPERATURE 
ONE     WALL    INSULATED 
X*  0.250  L  =0.0156 
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TEMPERATURE     AND  VELOCITY     PROFILES 
ONE    WALL    CONSTANT     TEMPERATURE 
ONE    WALL   INSULATED 
X=l.00  L  =0.0625 
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FIGURE      230 
DIMENSIONLESS    HEAT     FLUX 
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FIGURE      25 

TEMPERATURE    AND  VELOCITY     PROFILES 
CONSTANT    WALL   TEMPERATURE 
(ASYMMETRIC    CASE) 

X=  0.005  L  =0.0003125 
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FIGURE    26 

TEMPERATURE     AND   VELOCITY      PROFILES 
CONSTANT    WALL    TEMPERATURE 
(ASYMMETRIC   CASE) 
X=  0.050  L  =0.003125 
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FIGURE    27 

TEMPERATURE    AND    VELOCITY    PROFILES 
CONSTANT    WALL    TEMPERATURE 

.      (ASYMMETRIC  CASE) 

X=  0.250  L    =0.0156 
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FIGURE     28 

TEMPERATURE    AND   VELOCITY      PROFILES 
CONSTANT     WALL     TEMPERATURE 
(ASYMMETRIC   CASE) 

X*  1.00  L= 0.0625 
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CONSTANT    WALL    TEMPERATURE 
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..JOB0191F,LUNDBERG   ' 
PROGRAM  HEATX23. 
C   ,   CONSTANT  HEAT  INPUT 
C      BODOIA  VELOCITY  DATA 

C      THIS  PROGRAM  USES   A  REDUCED  GRID   OF.  X= .000 1 ♦ Y  =  . 05 

C      FOR  THE  FIRST  20  STREAMWISE  CALCULATIONS. 

C      GRAPH  OUTPUTS  OF  NUSSELT  NUMBER  VERSUS  LENGTH  PARAMETER 

C      FOR'PRANDTL  NUMBERS  OF  0.  5  ,0.7 » 1.0 , 1 .6, 3. 2 ,AND  10.  '• 

C      ARE  OBTAINED  TO  50  PER  CENT  OF  THE  HYDRODYNAMIC 

C      DEVELOPMENT.   PRINT  OUTPUT  INCLUDES  NUSSELT  NUMBER* 

C      SOPE  OF  THE  TEMPERATURE  PROFILE  AT  THE  WALL  ,  AND  THE 

C  .    SLOPE  OF  THE  TEMPERATURE  PROFILE  AT  THE  WALL, AND  THE  MEAN 

C      MIXED  MEAN  TEMPERATURE.   IN  ADDITION,  VELOCITY  AND 

C      TEMPERATURE  VALUES  ARE  PRINTED  AT  EACH  GRID  POINT  AT 

.C '     SELECTED  STREAMWISE  LOCATIONS. 

DIMENSION  TOO, 23)  ,U(30, 23), V(30,23)  ,Y<23)  ,X(500)   .. 
1,XN(500) ,ITITLE(12) 

READ  500, ( ITITLE< I ),  1=1,6) 

READ  500,{ IT!TLE< I ) ,  1=7,12) 
500  FORMAT  (6A8) 

DO  750  KK=1,6 

GO  TO(40,41,42»43,44,45) »KK 

40  PR=0.7 

GO  TO  46   / 

41  PR=0.5 

GO  TO  888 

GO  TO  46  •  • 

42  PR=1.0 
GO  TO  46 

■  43  PR=1.6 

GO  TO  46 

44  PR=3.2 
GO  TO  46 

45  PR=10. 

46  PRINT  523 
523  FORMAT  (1H1) 

PRINT  17,  'PR 
17  FORMAT  (/////10X,17H  PRANDTL  NUMBER  =»F5.3) 
L=l 
GO  TO  150 

700  READ  701,  ((U(ItJ)t  J=13»22)»  1=2,22) 

701  FORMAT   (10F6.4) 
DO  710  1=2,22 

DO  711  J=l,12 
711  U(  I  ,J)=U(  1,1-3)  ' 
710  CONTINUE 
•  -•    DELX=.0001 

DELY=0.05 

R=1000. 

DO  110  L=2,21 
150  UTM=0.0 

UM=0.0  .  \         ,  •  • 

XL  =  L 

X(L)=(XL-1.0)*0.0001 

55 


1  !■'  ,-ft "7 TrT-riT-,--;Y<~ "^" **?. '  r7T-"7 5 1  '-nrrn* 


108 
109 


113 


112 

114 
111 


131 
118 


119 


120 


121 


DO  111  J=3»22 

XJ  =  J 

Y(.l)  =  <XJ-3.0)*0.05 

Y(23)=1.0 

U( 1,J)=1. 

U(l»23)=1.0 

V(1,J-1)=0.0 

V(1,J)=0. 

T(1,J)=1.0  *  *' 

T(1,J+1)=1.0 

T(1,J-1)=1.0 

T(L,1)=T(L,5) 

U(L,1)=U(L»5) 

U<L,23>=0.0 

U(L,2)=U(L,4) 

V(L»3)=0.0 

V(L,2)=V(L,4) 

THERMAL  BOUNDARY  CONDITIONS 

T(L,23)=T<L,22)+.05 

T(L,2)=T(L,4) 

IF(L-l)  109, 109*108 

V(L,J)=  V(L»J-1)  +  DELY/<2.*DELX*R)  *  ( U ( L-l ♦ J ) -U( L+l # J ) > 

A=  DELX/  (U(L.J)*PR*DELY**2) 

V(L,23)=0.0 

B=  <DELX*R*V(L,J) )  /  ( 2 .*DELY*U ( L, J ) ) 

T(L+1,J)=(A-B)*T(L,J+1 )  +  U.-2.*A)*T(L.J)  +  (A+B)*T(L,J-]J 

IF  (J-3)  113,113,112  -  -*.. 

UTM=UTM+,5*U(L,J)*T(L»J) 

UM=UM+.5*U(L»J) 

GO  TO  114 

UTM=UTM+U(L»J)*T(L»J) 

UM=UM+U(L,J) 

TM=UTM/UM 

CONTINUE 

Tl=    (  T(L,23)-T(L,22)  )  /  DELY 

T2=     T(L,23)  -  TM 

XN(L)=4.0*T1/T2  , 

PRINT  118,  X(L) ,XN(L) ,T2,T1,TM 


FORMAT  (// 
1F8.5,5X,8H 
PRINT119,  <T(L,J) 
PRINT119,  (T(L,J) 
FORMAT  {/9H  TEMP, 
PRINT120,  •  (U(L,J) 
PRINT120, 
FORMAT  (/ 
PRINT121, 
PRINT121, 
FORMAT  (/ 


4H  X  =,F6.4,5X,14H 
SLOPE  =,F8.5,5X,5H 
J=l,12) 
J=13,23) 
=,12F7.4) 
J=l»12) 
J=13,23) 

=,12F7.4) 
J=l,12) 
J=13,23) 
9H  Y  VELv  =,12F7,4) 


NUSSELT  NO. 
TM  *,F8.5) 


=,F8.5,5X,8H  TW-TM  =  , 


(U(L,J) 
9H  X  VEL. 
(V(L,J) 
(V(L,J) 


IF(L-l)  700,700,110 
110  CONTINUE 

CALL  DRAW  < 21 »X ,XN,MOD»0 ,LAB» I T ITLE, .04 , 10. ,0 ,0 ,0 >0 . 6 » 

DO  720  N=2,ll 
M=2*N-1 


6»1»LAST) 


56 


771 

720 

85 

70 


U(  1,N)=U(11,M) 
U(2,N)=U<21,M) 
T(2,N)=T(21,M) 
PRINT  771,  Mi  U<2,N) »T(2»N) 
FORMAT  (I10.2F10.-5) 
CONTINUE 
L  =  2 

DO  10  M=2,250 
XM  =  M 
Y=0.0 
DELY=0.1 
DELX=0.001 
X(M)=DELX*XM 
I=L+1 
N3N=M 
N  =  L+1 

IF(M-  99)  72,73,73 
READ  1  (U(  I,U)  ,U(  ItlO)  »U(  1,9) 
1 ,U( I ,5) *U( I ,4) ,U< I »3) »U( I .2)  ) 
1  FORMAT  (10F6.4) 
73  IF(M-99  )  50,51,52 
IF(M-149)  51,53,54 
IF'(M-199)  53,55,55 
UU=XM-  99. 

U(N,2)  =1.4388+.000740*UU 
U(N,3)  =1.4292+.000674*UU 
U(N,4)  =1.3993+.000492*UU 
U(N,5)  =1.3454+.000238*UU 
U(N,6)  =1.2628-.000034*UU 
U<N,7)  =1.1467-. 000262*UU 
U(N,8)  =  .9938-.000410*UU 
U(N,9)  =  .8022-.000452*UU 
U(N,10)=  .5720-.000388*UU 
U(N,11)=  .3042-.000230*UU 
GO  TO  50 


82 
72 


52 
54 
51 


,U( I ,8)  ,U(  1,7)  ,U( I .6) 


53 


55 


UU=XM-149. 

U(N,2)  =1.4758+.000290*UU 
=1.4629+.000266*UU 
=1.4239+.000194*UU 
=1.3573+.000092*UU 
=1.2611-. 000014*UU 
=1.1336-. 000104*UU 
=  .9733-.000160*UU 
=  .7796-.000176*UU 
•5526-.000150*UU  * 
•2926-.000092*UU 


U(N,3) 

U(N,4) 

U(N,5) 

U(N,6) 

U(N,7) 

U(N,8)  = 

U(N,9)  = 

U(N,10)= 

U(N,11)= 

60  TO  50 

UU=XM-199. 

U(.N,2)  =1.4903+.  0000121*UU 

U(N,3)  =1.4762+.0000110*UU 

U(N,4)  =1.4336+.0000080*UU 

U(N,5)  =1.3619+. 0000026*UU 

U(N,6)  =1.2604-. 0000005*UU 

U(N,7)  =1.1284-. 0000042*UU 
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T" 


U(N,3)  =  .9653-.0000066*UU  . 

U(N,9)  =  .7708-.0000072*UU 

U(.N,10)  =  .5451-.0000063*UU 

U(N,11)-  .2880-.0000475*UU 
50  UTM=0.0 

UM=0,0 

DO  11  J=2,ll  *  . 

XJ  =  J 

YU)-DELY*XJ-.2 

U(L,12)=0.0 

V(L,2>=0.0 
•  V(L,1)=V( 1,3)  * 

•   Y(12)=1.0 

T(L,12).=  T(L,11)+0.1 

T(L,1)=T(L,3) 

U(L,1)=U(L,3) 
8  V(L,J)=  V(L,J-1)  +  DELY/(2.*DELX*R)  *  ( U( L-l , J )-U ( L+l , J ) ) 
>=  DELX/  (U(LiJ)*PR*DELY**2) 

V(L,12)=0.0  , 

B=  (DELX*R*V(L,J) )  /  ( 2 .*DELY*U ( L, J ) > 
.   T(L+1,J)=(A-B)*T(L,J+1)+(1.~2.*A)*T(L»J)+(A+B)*T(L,J-1) 

IFU-2)  13,13,12 

13  UTM=UTM+.5*U(L»J)*T(L»J) 
UM=UM+.5*U(L»J) 

GO  TO  14 
12.  UTM=UTM+U(L,J)*T(L,J)  •  ,   > 

UM=UM+U(L,J)  X 

14  TM=UTM/UM 
11  CONTINUE 

Tl=   (  T(L,12)-T(L,11) )  /  DELY 
T2=    T(L,12)-TM 
XN(M)=4.0*T1/T2 
IF(N3N-50  )  31,31,30 

30  IF( (N3N/25)*25-N3N)  15,31,31 

31  PRINT  18,X(M) ,XN(M) ,T2,T1,TM 

18  FORMAT  (//  4H  X  = »F6. 3 »5X ♦ 14H  NUSSELT  NO.  =,F8.5»5X»8H  TW-TM  =, 
'  1F8.5,5X,8H  SLOPE  =,F8.5,5X,5H  TM  *»F8.5) 

PRINT  19,  (T(L,J) ♦  J=l,12) 

19  FORMAT  (/9H  TEMP.   =,12F7.4) 
PRINT  20,  (U(L,J) ♦  J=l,12) 

20  FORMAT  (/  9H  X  VEL.  =,12F7.4) 
PRINT  21,  (V(L,J)»  J=l,12) 

21  FORMAT  (/  9H  Y  VEL.  -,12F7.4) 

15  DO  60  J=2,ll 
T(L,J)=T(L+1»J) 
U(L-1,J)=U(L,J) 

60  U(L,J)=U(L+1»J) 
10  CONTINUE        '■ 

IF(KK-6)  760,761,761 

761  MOD=3 

GO  TO  762 
760  MOD=2 

762  CALL  DRAW(497»X,XN,MOD»0,LAB»ITITLE,.04,10.,0,0,0»0,6»  6»  It  LAST). 
750  CONTINUE 
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